# ==============================================================================
# ISEA4H离散全球网格系统（DGGS）涡旋检测脚本
# 核心功能：基于ADT（绝对动力地形）数据，通过ISEA4H网格转换、极值点检测、
#          区域扩展、属性计算等步骤，实现海洋涡旋（气旋/反气旋）的自动检测与保存
# ==============================================================================

# --------------------------
# 1. 依赖库导入（标注库功能用途）
# --------------------------
import h5netcdf  # 读写NetCDF格式数据文件
import numpy as np  # 数值计算（数组处理、数学运算）
import numba  # 即时编译（JIT）加速计算密集型函数
import logging  # 日志记录（调试、运行信息输出）
from collections import deque  # 双端队列（用于BFS算法扩展涡旋区域）
from scipy.ndimage import gaussian_filter1d  # 高斯平滑（用于涡旋轮廓平滑处理）
import os  # 文件/目录路径操作
import re  # 正则表达式（用于从文件名提取日期）
import PyDGGS  # ISEA4H离散全球网格系统核心库（网格编码/解码、邻居查询）

# --------------------------
# 2. 日志配置（统一日志格式与级别）
# --------------------------
logging.basicConfig(
    level=logging.INFO,  # 日志级别：INFO（输出关键运行信息）
    format='%(asctime)s - %(name)s - %(levelname)s - %(message)s',  # 日志格式：时间-模块名-级别-内容
    datefmt='%Y-%m-%d %H:%M:%S'  # 时间格式：年-月-日 时:分:秒
)
logger = logging.getLogger('ISEA4HEddyDetector')  # 日志实例命名（区分不同模块）

# --------------------------
# 3. 核心参数配置（标注参数含义与作用）
# --------------------------
dggs_level = 8  # ISEA4H网格分辨率级别（级别越高，网格单元越小，分辨率越高）
grid_type = PyDGGS.DGGSGridType_DGGS_ISEA4H  # 网格类型：ISEA4H离散全球网格
element_type = PyDGGS.DGGSElementType_Cell  # 网格元素类型：单元格（Cell）
topo_type = PyDGGS.DGGSTopoType_RelativeCells  # 拓扑查询类型：相对单元格（获取邻居）
kernel = PyDGGS.CreateKernelObj(grid_type)  # 创建ISEA4H网格核心计算实例
CELL_RADIUS_KM = 17.3  # 单个ISEA4H网格单元的半径（单位：km，对应level=8的分辨率）
MIN_EDDY_RADIUS_KM = 50  # 涡旋最小有效半径（过滤过小的无效涡旋，单位：km）
MAX_EDDY_RADIUS_KM = 300  # 涡旋最大有效半径（过滤过大的异常涡旋，单位：km）
MAX_EXPAND_STEPS = 8  # 涡旋区域BFS扩展的最大圈数（控制涡旋扩展范围）


# --------------------------
# 4. 辅助函数（标注输入输出与核心逻辑）
# --------------------------
def get_hash_key(code):
    """
    将PyDGGS的DGGSElementCode转换为可哈希的元组（用于字典键/集合元素）
    PyDGGS原生DGGSElementCode对象不可哈希，需提取核心属性（行、列、基础分区）构建元组
    
    参数：
        code: PyDGGS.DGGSElementCode - 网格单元的编码对象
    返回：
        tuple: (code.row, code.col, code.basePartition) - 可哈希的单元标识
    """
    return (code.row, code.col, code.basePartition)


# --------------------------
# 5. 涡旋属性类（标注每个属性的物理意义）
# --------------------------
class Eddy:
    """
    涡旋属性容器类：存储单个涡旋的所有物理特征与几何属性
    每个属性对应涡旋的一项关键指标，便于后续分析与保存
    """

    def __init__(self):
        # 1. 涡旋核心物理属性
        self.amplitude = 0  # 涡旋振幅（ADT极值与周围均值的差值，单位：原始ADT缩放前）
        self.effective_area = 0.0  # 有效面积（涡旋覆盖的实际面积，单位：m²）
        self.effective_contour_height = 0.0  # 有效轮廓ADT高度（轮廓处的ADT值，单位：m）
        self.effective_radius = 0  # 有效半径（涡旋中心到轮廓的平均距离，单位：km）
        self.inner_contour_height = 0.0  # 内部轮廓ADT高度（涡旋中心附近的ADT值，单位：m）
        
        # 2. 涡旋几何位置属性
        self.latitude = 0.0  # 涡旋中心纬度（单位：度）
        self.latitude_max = 0.0  # 涡旋极值点纬度（ADT最高/最低点的纬度，单位：度）
        self.longitude = 0.0  # 涡旋中心经度（单位：度）
        self.longitude_max = 0.0  # 涡旋极值点经度（ADT最高/最低点的经度，单位：度）
        
        # 3. 涡旋轮廓属性（存储为int16以节省空间，实际值需除以100）
        self.effective_contour_latitude = np.zeros(50, dtype=np.int16)  # 有效轮廓纬度（×100，单位：度）
        self.effective_contour_longitude = np.zeros(50, dtype=np.int16)  # 有效轮廓经度（×100，单位：度）
        self.effective_contour_shape_error = 0  # 有效轮廓形状误差（描述轮廓与圆形的偏差）
        self.speed_contour_latitude = np.zeros(50, dtype=np.int16)  # 速度轮廓纬度（×100，单位：度）
        self.speed_contour_longitude = np.zeros(50, dtype=np.int16)  # 速度轮廓经度（×100，单位：度）
        self.speed_contour_shape_error = 0  # 速度轮廓形状误差
        
        # 4. 涡旋结构与统计属性
        self.num_contours = 0  # 轮廓数量（涡旋的嵌套轮廓层数）
        self.num_point_e = 0  # 有效点数量（涡旋内有效网格单元数）
        self.num_point_s = 0  # 速度点数量（与有效点数量一致，兼容输出格式）
        self.speed_area = 0.0  # 速度面积（基于速度轮廓的面积，单位：m²）
        self.speed_average = 0  # 平均速度（涡旋内的平均流场速度，单位：原始缩放前）
        self.speed_contour_height = 0.0  # 速度轮廓ADT高度（速度轮廓处的ADT值，单位：m）
        self.speed_radius = 0  # 速度半径（基于速度轮廓的半径，单位：km）
        
        # 5. 时间与类型属性
        self.time = 0  # 涡旋检测时间（对应原始数据的时间戳）
        self.eddy_type = ''  # 涡旋类型（'high'=反气旋，'low'=气旋）
        
        # 6. 网格单元属性
        self.cells = []  # 涡旋包含的ISEA4H网格单元（存储单元的哈希键）
        self.uavg_profile = np.zeros(50, dtype=np.uint16)  # 平均速度剖面（随半径变化的速度分布）


# --------------------------
# 6. 加速计算函数（标注JIT加速与极值检测逻辑）
# --------------------------
@numba.jit(nopython=True)
def find_adt_extrema(adt_array):
    """
    （Numba JIT加速）查找ADT数据的局部极值点（高值=反气旋中心，低值=气旋中心）
    采用3×3窗口对比：中心值大于周围所有值则为高值极值，小于则为低值极值
    
    参数：
        adt_array: np.ndarray - 2D ADT数据数组（纬度×经度，含NaN填充的无效区域）
    返回：
        list: 极值点列表，每个元素为(i, j, type)，其中：
            i: 纬度索引，j: 经度索引，type: 'high'（高值）/'low'（低值）
    """
    extrema = []
    nlat, nlon = adt_array.shape  # 获取ADT数组的纬度、经度维度

    # 遍历除边界外的所有网格点（边界无完整3×3窗口，跳过）
    for i in range(1, nlat - 1):
        for j in range(1, nlon - 1):
            val = adt_array[i, j]
            if np.isnan(val):  # 跳过NaN（无效区域）
                continue

            # 提取3×3窗口，将中心值设为NaN（避免与自身对比）
            window = adt_array[i - 1:i + 2, j - 1:j + 2].copy()
            window[1, 1] = np.nan

            # 判断极值类型：大于窗口最大值→高值（反气旋），小于窗口最小值→低值（气旋）
            if val > np.nanmax(window):
                extrema.append((i, j, 'high'))
            elif val < np.nanmin(window):
                extrema.append((i, j, 'low'))

    return extrema


# --------------------------
# 7. 核心涡旋检测类（分方法标注功能与流程）
# --------------------------
class ISEA4HEddyDetector:
    """
    ISEA4H网格涡旋检测器：整合数据加载、网格转换、极值检测、区域扩展、
    属性计算、结果过滤与保存的完整涡旋检测流程
    """

    def __init__(self, resolution=8):
        """
        初始化涡旋检测器：设置核心参数，默认值基于ISEA4H网格特性与涡旋检测经验
        
        参数：
            resolution: int - ISEA4H网格分辨率级别（默认8，对应CELL_RADIUS_KM=17.3km）
        """
        # 1. 网格与数据基础参数
        self.resolution = resolution  # ISEA4H网格级别
        self.fill_value = -2147483647  # ADT数据中的无效值填充（NC文件标准填充值）
        self.adt_scale = 1e-5  # ADT数据默认缩放因子（多数NC文件中ADT以1e-5为单位存储）
        self.earth_radius = 6371e3  # 地球半径（单位：m，用于计算球面距离）
        
        # 2. 涡旋过滤参数
        self.min_eddy_size = 8  # 涡旋最小包含单元数（过滤过小的噪声区域）
        self.max_expand_steps = MAX_EXPAND_STEPS  # BFS最大扩展圈数
        self.min_amplitude = 0.05  # 涡旋最小振幅（过滤弱信号涡旋，单位：m）
        self.min_radius = MIN_EDDY_RADIUS_KM  # 涡旋最小半径（单位：km）
        self.max_radius = MAX_EDDY_RADIUS_KM  # 涡旋最大半径（单位：km）
        
        # 3. 轮廓处理参数
        self.contour_sigma = 1.0  # 高斯平滑系数（值越小保留轮廓细节越多，越大越平滑）
        
        # 4. 统计信息（记录检测过程中的关键数量）
        self.stats = {
            'initial_extrema': 0,  # 初始检测到的极值点总数
            'filtered_extrema': 0,  # 过滤掉的无效极值点数量
            'merged_eddies': 0,     # 合并的重叠涡旋数量
            'filtered_eddies': 0,   # 过滤掉的无效涡旋数量
            'final_eddies': 0,      # 最终保留的涡旋数量
            'cyclonic_eddies': 0,   # 气旋涡旋数量（low型）
            'anticyclonic_eddies': 0# 反气旋涡旋数量（high型）
        }

    def get_neighbors(self, cell_key):
        """
        获取指定ISEA4H网格单元的所有邻居单元（基于PyDGGS拓扑查询）
        
        参数：
            cell_key: tuple - 网格单元的哈希键（由get_hash_key生成）
        返回：
            PyDGGS.QueryResult - 邻居单元的查询结果（含所有邻居的DGGSElementCode）
        """
        # 1. 从哈希键重建DGGSElementCode对象
        code = PyDGGS.DGGSElementCode()
        code.row, code.col, code.basePartition = cell_key  # 恢复行、列、基础分区
        code.level = dggs_level  # 设置网格级别
        code.elementType = element_type  # 设置元素类型（单元格）
        
        # 2. 调用PyDGGS接口查询邻居
        query_result = PyDGGS.QueryResult()
        kernel.query(topo_type, code, query_result)  # 拓扑查询：获取相对邻居
        
        return query_result

    def load_data(self, nc_path):
        """
        加载NetCDF文件中的ADT数据与元信息（时间、经纬度），并预处理无效值
        
        参数：
            nc_path: str - NetCDF文件的完整路径
        """
        logger.info(f"开始加载数据: {nc_path}")
        with h5netcdf.File(nc_path, 'r') as f:
            # 1. 读取ADT数据并处理无效值
            adt_raw = f['adt'][:].squeeze()  # 读取ADT数据，挤压单维度（如时间维度）
            # 将填充值替换为NaN，并用缩放因子转换为实际ADT值（单位：m）
            self.adt = np.where(adt_raw == self.fill_value, np.nan, adt_raw) * self.adt_scale
            
            # 2. 读取经纬度与时间元信息
            self.lats = f['latitude'][:].astype(np.float64)  # 纬度数组（单位：度）
            self.lons = f['longitude'][:].astype(np.float64)  # 经度数组（单位：度）
            self.time = f['time'][0].astype(np.uint32)  # 时间戳（取第一个时间点，假设单时次数据）
            
            # 3. 从NC文件属性更新ADT缩放因子（若文件中定义）
            if 'scale_factor' in f['adt'].attrs:
                self.adt_scale = f['adt'].attrs['scale_factor']
                logger.info(f"从文件获取ADT缩放因子: {self.adt_scale}")

        # 4. 从文件名提取检测日期（用于结果文件名）
        self.file_date = self._extract_date_from_filename(nc_path)
        logger.info(f"从文件名提取的日期: {self.file_date}")
        logger.info(f"数据加载完成: 纬度点数={len(self.lats)}, 经度点数={len(self.lons)}")

    def _extract_date_from_filename(self, filename):
        """
        从文件名中提取日期（匹配YYYYMMDD格式，如"adt_20231001.nc"提取为"2023-10-01"）
        
        参数：
            filename: str - 文件名（含路径或不含路径）
        返回：
            str - 格式化日期（YYYY-MM-DD），若匹配失败则返回当前日期
        """
        # 提取纯文件名（去掉路径），用正则匹配8位数字（YYYYMMDD）
        match = re.search(r'(\d{4})(\d{2})(\d{2})', os.path.basename(filename))
        if match:
            # 分组匹配：第1组=年，第2组=月，第3组=日，格式化为YYYY-MM-DD
            return f"{match.group(1)}-{match.group(2)}-{match.group(3)}"
        else:
            # 匹配失败时返回当前日期（避免无日期标识）
            from datetime import datetime
            return datetime.now().strftime("%Y-%m-%d")

    def _grid_to_isea4h(self):
        """
        将经纬度网格的ADT有效点转换为ISEA4H网格单元，并建立单元与ADT值的映射（字典）
        作用：将传统经纬度数据转换为离散全球网格数据，便于后续涡旋区域扩展与拓扑分析
        """
        logger.info("开始转换经纬度网格为ISEA4H单元...")
        # 1. 生成经纬度网格（纬度×经度）
        lon_grid, lat_grid = np.meshgrid(self.lons, self.lats)
        # 展平为1D数组（便于遍历每个网格点）
        lon_flat = lon_grid.flatten()
        lat_flat = lat_grid.flatten()
        adt_flat = self.adt.flatten()

        # 2. 筛选有效ADT点（排除NaN）
        valid_mask = ~np.isnan(adt_flat)
        self.cell_adt = {}  # 键：ISEA4H单元哈希键，值：对应ADT值（单位：m）

        # 3. 遍历每个有效点，编码为ISEA4H单元并存储
        for lon, lat, adt in zip(lon_flat[valid_mask], lat_flat[valid_mask], adt_flat[valid_mask]):
            # 构建地理坐标对象
            coord = PyDGGS.GeoCoord()
            coord.lon = lon  # 经度（度）
            coord.lat = lat  # 纬度（度）
            # 编码为ISEA4H单元
            cellCode = PyDGGS.DGGSElementCode()
            cellCode.level = dggs_level  # 网格级别
            cellCode.elementType = element_type  # 元素类型（单元格）
            kernel.encode(coord, cellCode)  # 经纬度→ISEA4H单元编码
            # 转换为可哈希键，存入字典
            cell_id = get_hash_key(cellCode)
            self.cell_adt[cell_id] = adt

        logger.info(f"转换完成: 有效网格点={len(self.cell_adt)}")

    def detect_eddies(self):
        """
        完整涡旋检测流程：极值点检测→区域扩展→属性计算→合并重叠涡旋→过滤无效涡旋
        
        返回：
            list[Eddy] - 最终检测到的有效涡旋列表
        """
        logger.info("开始涡旋检测流程...")
        # 步骤1：将经纬度网格转换为ISEA4H单元
        self._grid_to_isea4h()

        # 步骤2：查找ADT数据中的局部极值点（涡旋中心候选）
        logger.info("查找ADT极值点...")
        extrema = find_adt_extrema(self.adt)
        self.stats['initial_extrema'] = len(extrema)
        logger.info(f"找到{len(extrema)}个极值点")

        # 步骤3：基于极值点扩展涡旋区域，生成初始涡旋
        logger.info("开始扩展涡旋区域...")
        eddies = []  # 初始涡旋列表
        processed_cells = set()  # 已处理的网格单元（避免重复扩展）
        filtered_count = 0  # 过滤的无效极值点数量

        # 遍历每个极值点，尝试扩展为涡旋
        for i, (i_ext, j_ext, eddy_type) in enumerate(extrema):
            # 每处理1000个极值点输出进度（避免长时间无反馈）
            if i % 1000 == 0 and i > 0:
                logger.info(f"已处理{i}/{len(extrema)}个极值点，当前找到{len(eddies)}个涡旋")

            # 提取极值点的经纬度（涡旋中心候选位置）
            center_lat = self.lats[i_ext]
            center_lon = self.lons[j_ext]

            # 过滤经度超出合理范围的极值点（-180~180度）
            if not (-180 <= center_lon <= 180):
                filtered_count += 1
                continue

            # 将极值点经纬度编码为ISEA4H单元（涡旋中心单元）
            coord = PyDGGS.GeoCoord()
            coord.lon = center_lon
            coord.lat = center_lat
            center_cell = PyDGGS.DGGSElementCode()
            center_cell.level = dggs_level
            center_cell.elementType = element_type
            kernel.encode(coord, center_cell)
            center_index = get_hash_key(center_cell)

            # 过滤已处理或无效的中心单元（避免重复或无ADT值的单元）
            if center_index in processed_cells or center_index not in self.cell_adt:
                filtered_count += 1
                continue

            # BFS扩展涡旋区域（从中心单元向外扩展）
            eddy_cells = self._expand_eddy(center_index, eddy_type, processed_cells)
            # 过滤单元数过少的无效涡旋（小于最小单元数）
            if len(eddy_cells) < self.min_eddy_size:
                filtered_count += 1
                continue

            # 计算涡旋的完整属性（半径、面积、轮廓等）
            eddy = self._calculate_eddy_properties(eddy_cells, center_lon, center_lat, eddy_type)
            if eddy:  # 计算成功则加入初始涡旋列表
                eddies.append(eddy)

        # 更新统计信息：过滤的极值点数量
        self.stats['filtered_extrema'] = filtered_count
        logger.info(f"涡旋扩展完成: 初始检测到{len(eddies)}个涡旋，过滤掉{filtered_count}个极值点")

        # 步骤4：合并重叠的涡旋（同一类型且单元重叠的涡旋合并为一个）
        # 注：当前代码中合并逻辑已注释，若需启用可取消注释
        # logger.info("开始合并相交涡旋...")
        # initial_eddy_count = len(eddies)
        # merged_eddies = self._merge_overlapping_eddies(eddies)
        # merged_count = initial_eddy_count - len(merged_eddies)
        # self.stats['merged_eddies'] = merged_count
        # logger.info(f"涡旋合并完成: 合并了{merged_count}个，剩余{len(merged_eddies)}个")
        merged_eddies = eddies  # 暂用初始涡旋列表（未合并）

        # 步骤5：过滤不合理的涡旋（如经度超出范围）
        logger.info("开始过滤不合理涡旋...")
        filtered_eddies = self._filter_eddies(merged_eddies)
        filtered_count = len(merged_eddies) - len(filtered_eddies)
        self.stats['filtered_eddies'] = filtered_count
        logger.info(f"涡旋过滤完成: 过滤掉{filtered_count}个，剩余{len(filtered_eddies)}个")

        # 步骤6：统计涡旋类型（气旋/反气旋）与最终数量
        self.stats['cyclonic_eddies'] = sum(1 for e in filtered_eddies if e.eddy_type == 'low')
        self.stats['anticyclonic_eddies'] = sum(1 for e in filtered_eddies if e.eddy_type == 'high')
        self.stats['final_eddies'] = len(filtered_eddies)

        # 输出检测统计摘要（便于用户快速了解结果）
        logger.info("===== 涡旋检测统计摘要 =====")
        logger.info(f"初始极值点: {self.stats['initial_extrema']}")
        logger.info(f"过滤的极值点: {self.stats['filtered_extrema']}")
        logger.info(f"合并的涡旋: {self.stats['merged_eddies']}")
        logger.info(f"过滤的涡旋: {self.stats['filtered_eddies']}")
        logger.info(f"最终检测到的涡旋: {self.stats['final_eddies']}")
        logger.info(f"气旋涡旋: {self.stats['cyclonic_eddies']}")
        logger.info(f"反气旋涡旋: {self.stats['anticyclonic_eddies']}")
        logger.info("==========================")

        return filtered_eddies

    def _expand_eddy(self, center_cell, eddy_type, processed_cells):
        """
        基于BFS（广度优先搜索）从中心单元向外扩展涡旋区域，控制扩展范围与条件
        
        参数：
            center_cell: tuple - 涡旋中心单元的哈希键
            eddy_type: str - 涡旋类型（'high'=反气旋，'low'=气旋）
            processed_cells: set - 已处理的单元集合（避免重复处理）
        返回：
            list[tuple] - 涡旋包含的所有网格单元（哈希键）
        """
        # 获取中心单元的ADT值（作为扩展的参考阈值）
        center_adt = self.cell_adt[center_cell]
        # 初始化BFS队列：从中心单元开始
        queue = deque([center_cell])
        eddy_cells = [center_cell]  # 存储涡旋包含的单元
        processed_cells.add(center_cell)  # 标记中心单元为已处理
        expand_counts = {center_cell: 0}  # 记录每个单元的扩展圈数（从中心开始的距离）

        # BFS循环：直到队列为空或涡旋单元数超出上限
        while queue and len(eddy_cells) < self.min_eddy_size * 10:
            current_cell = queue.popleft()  # 取出当前单元
            current_step = expand_counts[current_cell]  # 当前单元的扩展圈数

            # 超出最大扩展圈数则停止扩展（控制涡旋大小）
            if current_step >= self.max_expand_steps:
                continue

            # 获取当前单元的所有邻居
            neighbors = self.get_neighbors(current_cell)
            neighbor_count = neighbors.getSize()  # 邻居单元数量

            # 遍历每个邻居单元，判断是否加入涡旋
            for i in range(neighbor_count):
                # 提取邻居单元的编码并转换为哈希键
                neighbor_code = neighbors.getGridCodeByIndex(i)
                neighbor_key = get_hash_key(neighbor_code)

                # 跳过已处理或无ADT值的邻居单元
                if neighbor_key in processed_cells or neighbor_key not in self.cell_adt:
                    continue

                # 获取邻居单元的ADT值
                neighbor_adt = self.cell_adt[neighbor_key]

                # 扩展条件：根据涡旋类型判断ADT值是否符合（反气旋：邻居ADT≤中心ADT；气旋：邻居ADT≥中心ADT）
                # 后期扩展（超过一半圈数）放宽条件，允许小概率加入（避免涡旋过早停止）
                expand_condition = (eddy_type == 'high' and neighbor_adt <= center_adt * 1.05) or \
                                  (eddy_type == 'low' and neighbor_adt >= center_adt * 0.95)
                if expand_condition or (current_step > self.max_expand_steps // 2 and np.random.random() < 0.1):
                    # 符合条件：标记为已处理，加入涡旋列表与队列
                    processed_cells.add(neighbor_key)
                    eddy_cells.append(neighbor_key)
                    queue.append(neighbor_key)
                    expand_counts[neighbor_key] = current_step + 1  # 记录扩展圈数

        return eddy_cells

    def _extract_contour(self, eddy_cells):
        """
        从涡旋包含的单元中提取边界轮廓单元（边界单元：至少一个邻居不在涡旋内）
        
        参数：
            eddy_cells: list[tuple] - 涡旋包含的单元列表（哈希键）
        返回：
            list[tuple] - 涡旋边界轮廓单元列表（哈希键）
        """
        eddy_cell_set = set(eddy_cells)  # 转换为集合，加速邻居判断
        contour_cells = []  # 存储轮廓单元

        # 遍历每个涡旋单元，判断是否为边界单元
        for cell in eddy_cells:
            # 获取当前单元的所有邻居
            neighbors = self.get_neighbors(cell)
            neighbor_count = neighbors.getSize()
            is_boundary = False  # 是否为边界单元的标记

            # 检查是否有邻居不在涡旋内
            for i in range(neighbor_count):
                neighbor_key = get_hash_key(neighbors.getGridCodeByIndex(i))
                if neighbor_key not in eddy_cell_set:
                    is_boundary = True
                    break  # 只要有一个邻居不在涡旋内，就是边界单元

            if is_boundary:
                contour_cells.append(cell)

        return contour_cells

    def _sort_contour_points(self, lons, lats, center_lon, center_lat):
        """
        按涡旋中心到轮廓点的角度排序轮廓点（顺时针或逆时针），使轮廓有序
        
        参数：
            lons: list[float] - 轮廓点经度列表（单位：度）
            lats: list[float] - 轮廓点纬度列表（单位：度）
            center_lon: float - 涡旋中心经度（单位：度）
            center_lat: float - 涡旋中心纬度（单位：度）
        返回：
            tuple(list[float], list[float]) - 排序后的轮廓点经度、纬度列表
        """
        if len(lons) < 2:  # 点数过少，无需排序
            return lons, lats

        # 计算每个轮廓点相对于中心的经纬度差（弧度）
        dlon = np.radians(lons - center_lon)
        dlat = np.radians(lats - center_lat)
        # 计算角度（arctan2(dlat, dlon)：从x轴正方向逆时针到点的角度）
        angles = np.arctan2(dlat, dlon)

        # 按角度排序（从小到大，对应逆时针方向）
        sorted_indices = np.argsort(angles)
        sorted_lons = np.array(lons)[sorted_indices].tolist()
        sorted_lats = np.array(lats)[sorted_indices].tolist()

        return sorted_lons, sorted_lats

    def _smooth_contour(self, lons, lats):
        """
        对涡旋轮廓进行高斯平滑（减少噪声，使轮廓更连续），同时保持轮廓闭合性
        
        参数：
            lons: list[float] - 排序后的轮廓点经度列表（单位：度）
            lats: list[float] - 排序后的轮廓点纬度列表（单位：度）
        返回：
            tuple(list[float], list[float]) - 平滑后的轮廓点经度、纬度列表
        """
        if len(lons) < 4:  # 点数过少（<4），平滑效果差，直接返回原轮廓
            return lons, lats

        # 转换为numpy数组便于处理
        lons = np.array(lons)
        lats = np.array(lats)

        # 填充首尾点（复制前2个和后2个点）：避免平滑时边界效应（轮廓两端变形）
        lons_padded = np.concatenate([[lons[-2], lons[-1]], lons, [lons[0], lons[1]]])
        lats_padded = np.concatenate([[lats[-2], lats[-1]], lats, [lats[0], lats[1]]])

        # 高斯平滑：sigma=1.0（轻微平滑，保留轮廓主要形状），mode='nearest'（边界值延伸）
        lons_smoothed = gaussian_filter1d(lons_padded, sigma=self.contour_sigma, mode='nearest')
        lats_smoothed = gaussian_filter1d(lats_padded, sigma=self.contour_sigma, mode='nearest')

        # 去除填充部分，恢复原始轮廓长度
        lons_smoothed = lons_smoothed[2:-2].tolist()
        lats_smoothed = lats_smoothed[2:-2].tolist()

        # 确保平滑后的轮廓首尾闭合（首尾点经纬度差异小于0.5度则强制闭合）
        if abs(lons_smoothed[0] - lons_smoothed[-1]) > 0.5:
            lons_smoothed[-1] = lons_smoothed[0]
        if abs(lats_smoothed[0] - lats_smoothed[-1]) > 0.5:
            lats_smoothed[-1] = lats_smoothed[0]

        return lons_smoothed, lats_smoothed

    def _get_contour_lonlat(self, contour_cells, center_lon, center_lat):
        """
        整合轮廓处理流程：从轮廓单元获取经纬度→排序→平滑，得到最终轮廓点
        
        参数：
            contour_cells: list[tuple] - 轮廓单元列表（哈希键）
            center_lon: float - 涡旋中心经度（单位：度）
            center_lat: float - 涡旋中心纬度（单位：度）
        返回：
            tuple(list[float], list[float]) - 最终轮廓点经度、纬度列表
        """
        lons, lats = [], []
        # 遍历每个轮廓单元，解码为经纬度（取单元中心经纬度）
        for cell in contour_cells:
            # 从哈希键重建DGGSElementCode
            code = PyDGGS.DGGSElementCode()
            code.row, code.col, code.basePartition = cell
            code.level = dggs_level
            code.elementType = element_type

            # 解码单元为地理坐标（单元中心经纬度）
            coord = PyDGGS.GeoCoord()
            kernel.decode(code, coord)
            lons.append(coord.lon)
            lats.append(coord.lat)

        # 按角度排序轮廓点→高斯平滑，得到有序、平滑的轮廓
        sorted_lons, sorted_lats = self._sort_contour_points(lons, lats, center_lon, center_lat)
        smoothed_lons, smoothed_lats = self._smooth_contour(sorted_lons, sorted_lats)

        return smoothed_lons, smoothed_lats

    def _resample_contour(self, lons, lats, num_points=50):
        """
        对涡旋轮廓进行均匀采样，确保每个涡旋的轮廓点数量一致（默认50个），便于后续存储与分析
        
        参数：
            lons: list[float] - 平滑后的轮廓点经度列表（单位：度）
            lats: list[float] - 平滑后的轮廓点纬度列表（单位：度）
            num_points: int - 目标采样点数（默认50，与Eddy类中轮廓数组长度一致）
        返回：
            tuple(list[float], list[float]) - 均匀采样后的轮廓点经度、纬度列表
        """
        if not lons or len(lons) < 2:  # 无轮廓点或点数过少，返回全零数组
            return np.zeros(num_points).tolist(), np.zeros(num_points).tolist()

        # 计算轮廓相邻点之间的距离（单位：m，基于球面近似：1度≈111km）
        distances = []
        for i in range(len(lons) - 1):
            # 经度距离：考虑纬度校正（1度经度在不同纬度的实际距离不同）
            dx = (lons[i + 1] - lons[i]) * 111e3 * np.cos(np.radians(lats[i]))
            # 纬度距离：1度≈111km
            dy = (lats[i + 1] - lats[i]) * 111e3
            # 欧氏距离（近似球面距离）
            distances.append(np.sqrt(dx ** 2 + dy ** 2))

        total_length = sum(distances)  # 轮廓总长度（单位：m）
        if total_length <= 0:  # 轮廓长度为0（异常），返回全零数组
            return np.zeros(num_points).tolist(), np.zeros(num_points).tolist()

        # 生成等距离采样位置（从0到总长度，共num_points个点）
        sample_positions = np.linspace(0, total_length, num_points, endpoint=False)
        sampled_lons = []
        sampled_lats = []
        current_position = 0  # 当前遍历到的轮廓长度
        segment_idx = 0  # 当前遍历的轮廓段索引

        # 遍历每个采样位置，找到对应的轮廓点
        for pos in sample_positions:
            # 移动到包含当前采样位置的轮廓段
            while segment_idx < len(distances) and current_position + distances[segment_idx] < pos:
                current_position += distances[segment_idx]
                segment_idx += 1

            # 超出轮廓段范围（异常），用最后一个点填充
            if segment_idx >= len(distances):
                sampled_lons.append(lons[-1])
                sampled_lats.append(lats[-1])
                continue

            # 计算在当前段内的插值比例（0~1）
            segment_length = distances[segment_idx]
            frac = (pos - current_position) / segment_length if segment_length != 0 else 0.0

            # 线性插值得到采样点经纬度
            lon = lons[segment_idx] + frac * (lons[segment_idx + 1] - lons[segment_idx])
            lat = lats[segment_idx] + frac * (lats[segment_idx + 1] - lats[segment_idx])

            sampled_lons.append(lon)
            sampled_lats.append(lat)

        return sampled_lons, sampled_lats

    def _fill_contour_arrays(self, eddy, lon_contour, lat_contour, nb_sample=50):
        """
        将采样后的轮廓经纬度填充到Eddy对象的轮廓数组中（转换为int16以节省存储空间）
        
        参数：
            eddy: Eddy - 待填充的涡旋对象
            lon_contour: list[float] - 采样后的轮廓经度列表（单位：度）
            lat_contour: list[float] - 采样后的轮廓纬度列表（单位：度）
            nb_sample: int - 采样点数（默认50，需与Eddy类中轮廓数组长度一致）
        """
        # 均匀采样轮廓点（确保点数为nb_sample）
        sampled_lons, sampled_lats = self._resample_contour(lon_contour, lat_contour, nb_sample)

        # 转换为int16存储（×100以保留2位小数精度，减少存储空间）
        eddy.effective_contour_longitude[:] = np.array(sampled_lons, dtype=np.float32) * 100
        eddy.effective_contour_latitude[:] = np.array(sampled_lats, dtype=np.float32) * 100
        # 速度轮廓与有效轮廓一致（兼容输出格式）
        eddy.speed_contour_longitude = eddy.effective_contour_longitude.copy()
        eddy.speed_contour_latitude = eddy.effective_contour_latitude.copy()

    def _calculate_radius(self, eddy_cells):
        """
        计算涡旋的有效半径（涡旋中心到所有单元中心的平均距离，单位：km）
        
        参数：
            eddy_cells: list[tuple] - 涡旋包含的单元列表（哈希键）
        返回：
            int - 涡旋有效半径（单位：km，取整）
        """
        if not eddy_cells:  # 无单元，返回0
            return 0

        # 1. 获取所有涡旋单元的中心经纬度
        lons, lats = [], []
        for cell in eddy_cells:
            code = PyDGGS.DGGSElementCode()
            code.row, code.col, code.basePartition = cell
            code.level = dggs_level
            code.elementType = element_type

            coord = PyDGGS.GeoCoord()
            kernel.decode(code, coord)
            lons.append(coord.lon)
            lats.append(coord.lat)

        # 2. 计算涡旋中心（所有单元经纬度的均值）
        center_lon, center_lat = np.mean(lons), np.mean(lats)

        # 3. 计算每个单元到中心的球面距离（单位：km）
        distances = []
        for lon, lat in zip(lons, lats):
            # 转换为弧度
            dlat = np.radians(lat - center_lat)
            dlon = np.radians(lon - center_lon)
            # 半正矢公式计算球面距离（单位：m）
            a = np.sin(dlat / 2) ** 2 + np.cos(np.radians(center_lat)) * np.cos(np.radians(lat)) * np.sin(dlon / 2) ** 2
            c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
            distance_km = (self.earth_radius * c) / 1000  # 转换为km
            distances.append(distance_km)

        # 4. 平均距离作为有效半径（取整）
        return np.uint16(np.mean(distances)) if distances else 0

    def _calculate_area(self, eddy_cells):
        """
        基于涡旋有效半径计算涡旋面积（假设涡旋为圆形，面积=πr²，单位：m²）
        
        参数：
            eddy_cells: list[tuple] - 涡旋包含的单元列表（哈希键）
        返回：
            np.float32 - 涡旋有效面积（单位：m²）
        """
        radius_km = self._calculate_radius(eddy_cells)  # 获取有效半径（km）
        if radius_km == 0:  # 半径为0，面积为0
            return 0.0
        # 转换半径为米，计算圆面积（πr²）
        radius_m = radius_km * 1000
        return np.float32(np.pi * (radius_m) ** 2)

    def _calculate_velocity_profile(self, eddy_cells, center_lon, center_lat):
        """
        计算涡旋的平均速度剖面（随半径变化的速度分布，模拟实际流场速度趋势）
        
        参数：
            eddy_cells: list[tuple] - 涡旋包含的单元列表（哈希键）
            center_lon: float - 涡旋中心经度（单位：度）
            center_lat: float - 涡旋中心纬度（单位：度）
        返回：
            np.ndarray - 速度剖面数组（长度50， dtype=np.uint16，对应不同半径的速度）
        """
        # 1. 计算每个单元到中心的距离（单位：km）
        distances = []
        for cell in eddy_cells:
            code = PyDGGS.DGGSElementCode()
            code.row, code.col, code.basePartition = cell
            code.level = dggs_level
            code.elementType = element_type

            coord = PyDGGS.GeoCoord()
            kernel.decode(code, coord)
            lon, lat = coord.lon, coord.lat

            # 计算平面距离（近似球面距离，单位：km）
            dx = (lon - center_lon) * 111e3 * np.cos(np.radians(center_lat)) / 1000  # 转换为km
            dy = (lat - center_lat) * 111e3 / 1000  # 转换为km
            distances.append(np.sqrt(dx ** 2 + dy ** 2))

        if not distances:  # 无距离数据，返回全零剖面
            return np.zeros(50, dtype=np.uint16)

        # 2. 生成归一化半径（0~1，对应从中心到边缘）
        max_dist = max(distances)
        normalized_dist = np.linspace(0, 1, 50)  # 50个点，覆盖0到最大距离

        # 3. 模拟速度剖面（指数衰减：中心速度高，边缘低，符合实际涡旋流场）
        # 公式：v = exp(-(r²)*5) * 10000，确保值在0~65535（uint16范围）
        profile = np.exp(-(normalized_dist ** 2) * 5) * 10000
        # 裁剪到uint16范围（0~65535），转换为整数类型
        return np.clip(profile, 0, 65535).astype(np.uint16)

    def _calculate_eddy_properties(self, eddy_cells, center_lon, center_lat, eddy_type):
        """
        计算涡旋的完整物理与几何属性，填充到Eddy对象中
        
        参数：
            eddy_cells: list[tuple] - 涡旋包含的单元列表（哈希键）
            center_lon: float - 涡旋中心经度（单位：度）
            center_lat: float - 涡旋中心纬度（单位：度）
            eddy_type: str - 涡旋类型（'high'=反气旋，'low'=气旋）
        返回：
            Eddy - 填充好属性的涡旋对象（若计算失败则返回None）
        """
        eddy = Eddy()  # 初始化涡旋对象
        eddy.eddy_type = eddy_type  # 设置涡旋类型
        eddy.cells = eddy_cells  # 存储涡旋包含的单元
        eddy.latitude = center_lat  # 设置涡旋中心纬度
        eddy.longitude = center_lon  # 设置涡旋中心经度
        eddy.time = self.time  # 设置时间戳
        eddy.num_point_e = eddy.num_point_s = len(eddy_cells)  # 有效点数量=单元数

        # 1. 计算ADT相关属性
        adt_values = np.array([self.cell_adt[cid] for cid in eddy_cells])  # 所有单元的ADT值
        # 振幅：ADT最大值与最小值的差值（除以缩放因子，恢复原始编码值）
        eddy.amplitude = np.uint16(np.abs(np.max(adt_values) - np.min(adt_values)) / self.adt_scale)
        eddy.effective_contour_height = np.min(adt_values)  # 有效轮廓ADT高度（最小值）
        eddy.inner_contour_height = np.max(adt_values)  # 内部轮廓ADT高度（最大值）
        eddy.speed_contour_height = np.mean(adt_values)  # 速度轮廓ADT高度（平均值）

        # 2. 处理涡旋轮廓（提取→排序→平滑→采样→填充）
        contour_cells = self._extract_contour(eddy_cells)  # 提取轮廓单元
        lon_contour, lat_contour = self._get_contour_lonlat(contour_cells, center_lon, center_lat)  # 轮廓经纬度
        self._fill_contour_arrays(eddy, lon_contour, lat_contour, 50)  # 填充轮廓数组

        # 3. 计算半径与面积属性
        eddy.effective_radius = self._calculate_radius(eddy_cells)  # 有效半径（km）
        eddy.effective_area = self._calculate_area(eddy_cells)  # 有效面积（m²）
        eddy.speed_radius = int(eddy.effective_radius * 0.8)  # 速度半径（有效半径的80%）
        eddy.speed_area = eddy.effective_area * 0.64  # 速度面积（有效面积的64%，对应半径80%）

        # 4. 计算轮廓形状误差（简化计算：轮廓单元数/10，取整）
        eddy.effective_contour_shape_error = np.uint8(len(contour_cells) / 10)
        eddy.speed_contour_shape_error = np.uint8(eddy.effective_contour_shape_error * 0.9)  # 速度轮廓误差为90%

        # 5. 计算平均速度与速度剖面
        # 平均速度：ADT平均值/缩放因子×100（模拟速度编码值）
        eddy.speed_average = np.uint16(np.mean(adt_values) / self.adt_scale * 100)
        # 速度剖面：随半径变化的速度分布
        eddy.uavg_profile = self._calculate_velocity_profile(eddy_cells, center_lon, center_lat)

        # 6. 找到ADT极值点的经纬度（高值/低值对应的单元中心）
        if eddy_type == 'high':
            extrema_idx = np.argmax(adt_values)  # 高值索引（反气旋）
        else:
            extrema_idx = np.argmin(adt_values)  # 低值索引（气旋）
        extrema_cell = eddy_cells[extrema_idx]  # 极值点所在单元
        # 解码极值单元的经纬度
        code = PyDGGS.DGGSElementCode()
        code.row, code.col, code.basePartition = extrema_cell
        code.level = dggs_level
        coord = PyDGGS.GeoCoord()
        kernel.decode(code, coord)
        eddy.latitude_max = coord.lat  # 极值点纬度
        eddy.longitude_max = coord.lon  # 极值点经度

        # 7. 计算轮廓数量（简化：轮廓单元数//5，取整）
        eddy.num_contours = np.uint16(len(contour_cells) // 5)

        return eddy

    def _merge_overlapping_eddies(self, eddies):
        """
        合并重叠的涡旋（仅合并同一类型的涡旋，且单元有交集）
        
        参数：
            eddies: list[Eddy] - 待合并的涡旋列表
        返回：
            list[Eddy] - 合并后的涡旋列表
        """
        if not eddies:  # 无涡旋，直接返回空列表
            return []

        logger.info(f"开始合并{len(eddies)}个涡旋...")
        merged_eddies = []  # 存储合并后的涡旋
        merge_count = 0  # 合并的涡旋数量

        # 遍历每个待合并的涡旋
        for eddy in eddies:
            merged_idx = -1  # 待合并的目标涡旋索引（-1表示无匹配）
            # 检查是否与已合并的涡旋重叠且类型一致
            for j in range(len(merged_eddies)):
                m_eddy = merged_eddies[j]
                # 仅合并同一类型的涡旋（气旋只合并气旋，反气旋只合并反气旋）
                if m_eddy.eddy_type != eddy.eddy_type:
                    continue
                # 检查单元是否有交集（集合交集非空则重叠）
                if len(set(eddy.cells) & set(m_eddy.cells)) > 0:
                    merged_idx = j
                    break

            # 找到可合并的涡旋：合并到目标涡旋
            if merged_idx != -1:
                merged_eddies[merged_idx] = self._merge_two_eddies(merged_eddies[merged_idx], eddy)
                merge_count += 1
            else:
                # 无重叠涡旋：直接加入合并列表
                merged_eddies.append(eddy)

        logger.info(f"涡旋合并完成: 合并了{merge_count}个，剩余{len(merged_eddies)}个")
        return merged_eddies

    def _merge_two_eddies(self, e1, e2):
        """
        合并两个重叠的涡旋：保留振幅大的涡旋作为主涡旋，整合属性
        
        参数：
            e1: Eddy - 第一个待合并的涡旋
            e2: Eddy - 第二个待合并的涡旋
        返回：
            Eddy - 合并后的涡旋
        """
        # 选择振幅大的涡旋作为主涡旋（保留更强信号的涡旋属性）
        if e1.amplitude >= e2.amplitude:
            main_eddy, sub_eddy = e1, e2
        else:
            main_eddy, sub_eddy = e2, e1

        # 初始化合并后的涡旋
        merged = Eddy()
        merged.eddy_type = main_eddy.eddy_type  # 保留主涡旋类型
        merged.cells = list(set(main_eddy.cells + sub_eddy.cells))  # 合并单元（去重）
        merged.latitude = main_eddy.latitude  # 保留主涡旋中心纬度
        merged.longitude = main_eddy.longitude  # 保留主涡旋中心经度
        merged.amplitude = max(main_eddy.amplitude, sub_eddy.amplitude)  # 取最大振幅
        merged.effective_radius = max(main_eddy.effective_radius, sub_eddy.effective_radius)  # 取最大半径
        merged.time = main_eddy.time  # 保留主涡旋时间戳

        # 整合其他属性
        merged.num_point_e = merged.num_point_s = len(merged.cells)  # 单元数更新为合并后的数量
        merged.effective_area = main_eddy.effective_area + sub_eddy.effective_area * 0.5  # 面积：主面积+子面积的50%
        merged.speed_area = merged.effective_area * 0.64  # 速度面积：有效面积的64%

        return merged

    def _filter_eddies(self, eddies):
        """
        过滤无效涡旋：当前仅过滤经度超出-180~180度的涡旋（避免跨日期线的异常涡旋）
        
        参数：
            eddies: list[Eddy] - 待过滤的涡旋列表
        返回：
            list[Eddy] - 过滤后的有效涡旋列表
        """
        logger.info(f"开始过滤{len(eddies)}个涡旋...")
        filtered_eddies = []  # 存储有效涡旋
        filtered_count = 0  # 过滤的无效涡旋数量

        # 遍历每个涡旋，检查轮廓经度范围
        for eddy in eddies:
            # 提取轮廓经度（需除以100，恢复原始度值）
            contour_lons = eddy.effective_contour_longitude / 100
            lon_min = np.min(contour_lons)  # 轮廓最小经度
            lon_max = np.max(contour_lons)  # 轮廓最大经度

            # 过滤经度超出-180~180度的涡旋
            if lon_min < -180 or lon_max > 180:
                filtered_count += 1
                logger.debug(f"过滤涡旋: 经度范围 [{lon_min:.2f}, {lon_max:.2f}]")
                continue

            # 有效涡旋：加入列表
            filtered_eddies.append(eddy)

        logger.info(f"涡旋过滤完成: 保留{len(filtered_eddies)}个，过滤掉{filtered_count}个")
        return filtered_eddies

    def save_to_nc(self, eddies, output_path_base):
        """
        将检测到的涡旋按类型（气旋/反气旋）分别保存为NetCDF文件
        
        参数：
            eddies: list[Eddy] - 待保存的涡旋列表
            output_path_base: str - 输出文件路径前缀（不含类型与后缀）
        """
        if not eddies:  # 无涡旋，不生成文件
            logger.info("未检测到涡旋，不生成输出文件")
            return

        # 按类型分类涡旋（气旋=low，反气旋=high）
        cyclonic = [e for e in eddies if e.eddy_type == 'low']  # 气旋涡旋
        anticyclonic = [e for e in eddies if e.eddy_type == 'high']  # 反气旋涡旋

        # 分别保存两种类型的涡旋
        if anticyclonic:
            output_path = f"{output_path_base}_{self.file_date}_Anticyclonic.nc"
            self._save_eddy_group(anticyclonic, output_path)
        if cyclonic:
            output_path = f"{output_path_base}_{self.file_date}_Cyclonic.nc"
            self._save_eddy_group(cyclonic, output_path)

    def _save_eddy_group(self, eddies, output_path):
        """
        将单一类型的涡旋列表保存为NetCDF文件，定义变量与维度
        
        参数：
            eddies: list[Eddy] - 单一类型的涡旋列表（气旋或反气旋）
            output_path: str - 输出NetCDF文件的完整路径
        """
        logger.info(f"开始保存{len(eddies)}个涡旋到{output_path}")
        with h5netcdf.File(output_path, 'w') as f:
            # 定义NetCDF维度：obs=涡旋数量，NbSample=轮廓采样点数（50）
            f.dimensions['obs'] = len(eddies)
            f.dimensions['NbSample'] = 50

            # 定义变量列表：(变量名, 维度, 数据类型)
            variables = [
                ('amplitude', ('obs',), np.uint16),
                ('effective_area', ('obs',), np.float32),
                ('effective_contour_height', ('obs',), np.float32),
                ('effective_contour_latitude', ('obs', 'NbSample'), np.int16),
                ('effective_contour_longitude', ('obs', 'NbSample'), np.int16),
                ('effective_contour_shape_error', ('obs',), np.uint8),
                ('effective_radius', ('obs',), np.uint16),
                ('inner_contour_height', ('obs',), np.float32),
                ('latitude', ('obs',), np.float32),
                ('latitude_max', ('obs',), np.float32),
                ('longitude', ('obs',), np.float32),
                ('longitude_max', ('obs',), np.float32),
                ('num_contours', ('obs',), np.uint16),
                ('num_point_e', ('obs',), np.uint16),
                ('num_point_s', ('obs',), np.uint16),
                ('speed_area', ('obs',), np.float32),
                ('speed_average', ('obs',), np.uint16),
                ('speed_contour_height', ('obs',), np.float32),
                ('speed_contour_latitude', ('obs', 'NbSample'), np.int16),
                ('speed_contour_longitude', ('obs', 'NbSample'), np.int16),
                ('speed_contour_shape_error', ('obs',), np.uint8),
                ('speed_radius', ('obs',), np.uint16),
                ('time', ('obs',), np.uint32),
                ('uavg_profile', ('obs', 'NbSample'), np.uint16),
            ]

            # 创建变量并填充数据
            for name, dims, dtype in variables:
                # 创建NetCDF变量
                var = f.create_variable(name, dims, dtype)
                # 遍历每个涡旋，填充对应属性
                for i, eddy in enumerate(eddies):
                    var[i] = getattr(eddy, name)  # 从Eddy对象获取属性值

        logger.info(f"成功保存{len(eddies)}个涡旋到{output_path}")


# --------------------------
# 8. 批量处理函数（标注目录处理流程与日志配置）
# --------------------------
def process_directory(input_dir, output_dir, resolution=8):
    """
    批量处理指定目录下的所有NetCDF文件，检测涡旋并保存结果
    
    参数：
        input_dir: str - 输入目录（存放ADT NetCDF文件）
        output_dir: str - 输出目录（存放涡旋检测结果文件）
        resolution: int - ISEA4H网格分辨率级别（默认8）
    """
    # 创建输出目录（若不存在）
    os.makedirs(output_dir, exist_ok=True)

    # 添加文件日志处理器（记录批量处理日志到文件）
    file_handler = logging.FileHandler(os.path.join(output_dir, 'eddy_detection_batch.log'))
    file_handler.setFormatter(logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s'))
    logger.addHandler(file_handler)

    # 筛选目录下的NetCDF文件（仅后缀为.nc的文件）
    nc_files = [f for f in os.listdir(input_dir) if f.endswith('.nc')]
    logger.info(f"发现{len(nc_files)}个NetCDF文件")

    # 统计批量处理结果
    total_eddies = 0  # 总检测到的涡旋数量
    total_files = 0   # 成功处理的文件数量

    # 遍历每个NetCDF文件，执行检测
    for i, nc_file in enumerate(nc_files):
        try:
            # 输出当前处理进度
            logger.info(f"===== 处理文件 {i + 1}/{len(nc_files)}: {nc_file} =====")
            nc_path = os.path.join(input_dir, nc_file)  # 完整文件路径

            # 初始化涡旋检测器
            detector = ISEA4HEddyDetector(resolution=resolution)
            detector.min_amplitude = 0.07  # 调整最小振幅阈值（过滤弱涡旋）
            logger.info(f"开始检测涡旋，平滑系数={detector.contour_sigma}")

            # 执行检测流程：加载数据→检测涡旋
            detector.load_data(nc_path)
            eddies = detector.detect_eddies()

            # 输出当前文件的检测结果
            logger.info(f"===== 文件 {nc_file} 最终结果 =====")
            logger.info(f"检测到{len(eddies)}个涡旋")
            if eddies:  # 有涡旋时计算平均半径与振幅
                avg_radius = np.mean([e.effective_radius for e in eddies])
                avg_amplitude = np.mean([e.amplitude * detector.adt_scale for e in eddies])
                logger.info(f"平均半径: {avg_radius:.2f}km")
                logger.info(f"平均振幅: {avg_amplitude:.4f}m")
            
            # 保存当前文件的涡旋检测结果
            base_output_name = os.path.join(output_dir, os.path.splitext(nc_file)[0])
            detector.save_to_nc(eddies, base_output_name)

            # 更新统计信息
            total_eddies += len(eddies)
            total_files += 1

        except Exception as e:
            # 捕获并记录单个文件处理中的异常，不中断批量处理
            logger.error(f"处理文件 {nc_file} 时出错: {str(e)}", exc_info=True)

    # 批量处理完成后输出总结
    logger.info(f"===== 批量处理完成 =====")
    logger.info(f"共处理{total_files}个文件，检测到{total_eddies}个涡旋")
    logger.info("========================")

    # 移除文件日志处理器（避免影响后续日志）
    logger.removeHandler(file_handler)


# --------------------------
# 9. 程序入口（标注默认路径与执行方式）
# --------------------------
if __name__ == "__main__":
    # 默认输入输出路径（用户可根据实际数据路径修改）
    input_directory = "D:/EddyTrajData/py-eddy-tracker-master/tests/my_test/data"  # ADT数据文件目录
    output_directory = "./EddyExtract"  # 涡旋检测结果输出目录
    # 执行批量处理（分辨率级别=8）
    process_directory(input_directory, output_directory, resolution=8)